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Abstract 

Recent experimental studies investigating the neuronal regulation of rapid eye movement (REM) sleep have identified 
mutually inhibitory synaptic projections among REIVl sleep-promoting (REM-on) and REM sleep-inhibiting (REIVl-off) 
neuronal populations that act to maintain the REM sleep state and control its onset and offset. The control mechanism of 
mutually inhibitory synaptic interactions mirrors the proposed flip-flop switch for sleep-wake regulation consisting of 
mutually inhibitory synaptic projections between wake- and sleep-promoting neuronal populations. While a number of 
synaptic projections have been identified between these REM-on/REM-off populations and wake/sleep-promoting 
populations, the specific interactions that govern behavioral state transitions have not been completely determined. Using 
a minimal mathematical model, we investigated behavioral state transition dynamics dictated by a system of coupled flip- 
flops, one to control transitions between wake and sleep states, and another to control transitions into and out of REM 
sleep. The model describes the neurotransmitter-mediated inhibitory interactions between a wake- and sleep-promoting 
population, and between a REM-on and REM-off population. We proposed interactions between the wake/sleep and REM- 
on/REM-off flip-flops to replicate the behavioral state statistics and probabilities of behavioral state transitions measured 
from experimental recordings of rat sleep under ad libitum conditions and after 24 h of REM sleep deprivation. Reliable 
transitions from REM sleep to wake, as dictated by the data, indicated the necessity of an excitatory projection from the 
REM-on population to the wake-promoting population. To replicate the increase in REM-wake-REM transitions observed 
after 24 h REM sleep deprivation required that this excitatory projection promote transient activation of the wake- 
promoting population. Obtaining the reliable wake-nonREM sleep transitions observed in the data required that activity of 
the wake-promoting population modulated the interaction between the REM-on and REM-off populations. This analysis 
suggests neuronal processes to be targeted in further experimental studies of the regulatory mechanisms of REM sleep. 
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introduction 

Theories on the neuronal control for rapid eye movement 
(REM) sleep have been dominated by the cholinergic hypothesis 
(see [1] for review). Based on a wealth of experimental evidence 
collected since the identification of REM sleep in the 1950s [2] , 
the cholinergic hypothesis posits that the REM sleep state is 
initiated and maintained by the activity of cholinergic neurons in 
areas of the pons, including the laterodorsal and pedunculopontine 
tegmental nuclei (LDT/PPT). This hypothesis is synthesized in the 
reciprocal interaction model [3,4] for REM sleep in which regular 
transitions between REM and nonREM (NREM) sleep are 
generated by excitatory and inhibitory synaptic projections 
between the cholinergic REM-promoting (REM-on) LDT/PPT 
and monoaminergic, REM-suppressing (REM-ofI) neuronal pop- 
ulations including the locus coeruleus (EC) and the dorsal raphe 
(DR). 

Recent experimental results have challenged the cholinergic 
hypothesis for REM sleep control with the identification of several 
additional neuronal groups with REM-on and REM-off activity 



profiles [5-7]. In the rat, REM-on groups include the sublater- 
odorsal tegmental nucleus (SLD), portions of the ventrolaterEil 
periaqueductal gray matter (vlPAG), areas of the lateral hypothal- 
amus (EH) and the dorsal paragigantocellular nucleus (DPGi). 
Neuronal groups with REM-off activity profiles include other 
portions of the vlPAG and the dorsal part of the deep 
mesencephalic nucleus (dDPME). Based on the identification of 
these REM-associated neuronal groups and their synaptic 
interactions, alternative theories for the neuronal control of 
REM sleep have been proposed wherein GABA is the primary 
neurotransmitter and mutually inhibitory synaptic interactions 
govern activity of these groups and thus transitions of REM sleep 
[5,7]. However, debate continues as to which of these GABAergic 
populations are integrally responsible for generating transitions 
into and out of REM sleep. For example, Lu and colleagues [5] 
propose that REM regulation is controlled by a core REM-on/ 
REM-off flip-flop switch composed of mutually inhibitory 
(GABAergic) synaptic projections between REM-oflF neurons in 
the vlPAG and adjacent lateral pontine tegmentum (LPT), and 
REM-on neurons in the SLD. On the other hand, Luppi and 
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colleagues [8] propose that REM sleep transitions are controlled 
by a more distributed network of inhibitory projections among 
REM-ofT neurons in the vlPAG and dDPME, and REM-on 
neurons in the vlPAG, LH and DPGi. 

Another factor influencing REM sleep generation is the REM 
sleep homeostatic drive. Multiple studies have provided evidence 
that REM sleep is homeostaticaUy regulated independently from 
NREM sleep [9,10]. For example, increases in attempts to 
transition into REM sleep have been observed during REM sleep 
deprivation studies in which these transitions are prevented, and 
REM sleep rebound reliably occurs during recovery sleep from 
periods of REM sleep deprivation as short as 2 hours [10]. 
Additionally, anesthetics have differential effects in terms of 
satisfying the needs of NREM and REM sleep following sleep 
deprivation, suggesting separate mechanisms for NREM and 
REM sleep homeostasis. In particular, anesthesia induced by the 
inhaled anesthetic sevoflurane following total sleep deprivation 
eliminated homeostatic increases in NREM sleep but not in REM 
sleep [11]. The physiological substrate dictating REM sleep 
homeostasis has not been identified. Recent studies have suggested 
that melanin-concentrating hormone (MCH) [12] or the satiety 
molecule Nesfatin-1 [13] may be involved. Additionally, metabolic 
activation of neurons in preoptic areas of the hypothalamus has 
been found to be strongly related to homeostatic pressure for REM 
sleep in REM sleep deprivation studies [14]. Analysis of sleep 
patterns during REM rebound and recovery following sleep 
deprivation have suggested that REM sleep may be regulated on 
both short-term and long-term time scales [15,16]. The longer 
time-scale regulates the daily amount of REM sleep and the short- 
term process dictates transitions between NREM and REM sleep 
during sleep episodes [15]. 

Taken together, the proposed mutually inhibitory network of 
neuronal populations governing transitions of REM sleep, coupled 
with a REM sleep homeostatic drive, mirrors the proposed flip- 
flop switch for sleep regulation composed of mutually inhibitory 
synaptic projections between the wake-promoting monoaminergic 
populations locus coeruleus (LC) and dorsal raphe (DR) and the 
sleep-promoting GABAergic ventral lateral preoptic nucleus 
(VLPO), under the influence of the putative adenosine-mediated 
homeostatic: sleep drive [17]. How these two flip-flops may be 
coupled to produce transitions between wake, NREM and REM 
sleep observed in the majority of mammalian species is not 
completely determined. Synaptic projections from the wake- 
promoting LC and DR targeting REM-on and REM-off neuronal 
groups have been identified that act to suppress REM sleep 
[8,18,19]. However, limitations in recording simultaneously from 
these multiple areas at spatial and temporal resolutions to 
determine causal patterns of activity and interaction for transitions 
into and out of sleep and wake states hamper the ability to confirm 
or refute this and competing hypotheses of the sleep-wake 
regulatory network. 

Recently, physiologically-based mathematical models have been 
introduced as a means to test the different and competing 
hypotheses for the sleep-wake regulatory network [20-25]. In 
previous work, we used a modeling formalism for the neurotrans- 
mitter-mediated interactions among wake-, sleep- and REM-sleep 
neuronal populations to investigate the cholinergic hypothesis for 
REM sleep generation in rats [24]. Here, we apply the same 
modeling formalism and develop a coupled flip-flop model to 
investigate how a wake/sleep flip-flop and a REM-on/REM-off 
flip-flop can interact to produce accurate behavioral state 
transition dynamics for the rat. We use experimental recordings 
of rat sleep behavior under ad libitum (baseline) conditions and 
during REM sleep rebound after 24 h of REM sleep deprivation 



to motivate a minimal set of projections between the flip-flops to 
account for sleep-wake patterns in both conditions. 

The dynamics of a single mutually inhibitory flip-flop model 
have been weU studied in the context of governing the transitions 
between wake and sleep states under control of a homeostatic sleep 
drive [21,26,27]. As we and others have shown previously, the 
dynamics are those of movement around a hysteresis loop with the 
homeostatic drive increasing and decreasing through bifurcation 
points to induce transitions between low and high population 
activity levels. As we discuss below, hysteresis loop dynamics have 
an inherent symmetry and regularity that may seem incompatible 
with the highly variable sleep-wake activity patterns of the rat. As 
in our previous modeling study of rat sleep-wake patterns [24], we 
include in our coupled flip-flop model several physiologically 
motivated stochastic components that introduce noise in our 
model solutions. To understand the implications of underlying 
hysteresis loop dynamics, we present a detailed analysis of the 
effects of the stochastic components on simulated bout durations 
for a single flip-flop model. 

Results 

Dynamics of a single flip-flop 

In this section we briefly review the hysteresis loop dynamics of 
a single flip-flop and analyze the effects on these dynamics of 
physiologically motivated sources of noise. The goal of this analysis 
is to understand how the inherent symmetric and regular 
dynamics of a hysteresis loop can be modulated to simulate the 
variability of rat sleep-wake behavior, including the reported 
qualitative differences in the distributions of wake and sleep bout 
durations [28]. 

Hysteresis loop dynamics of a single flip-flop. We 

consider the wake/sleep flip-flop consisting of wake-promoting 
and sleep-promoting populations with reciprocal inhibitory 
neurotransmitter-mediated projections between them (Fig. lA). 
The dynamics of the REM-on/REM-oflf flip-flop are analogous. 
In the figure, rectangles represent neuronal groups and are labeled 
with their firing rate variables, fW{t) and fS(t), while circles 
represent neurotransmitter concentrations expressed by the 
neuronal groups and are labeled with their variable names, 
cW{t) and cS(t). The wake/sleep flip-flop represents the mutually 
inhibitory synaptic interactions between the wake -promoting locus 
coeruleus (LC), dorsal raphe (DR) and the tuberomammilary 
nucleus (TMN) (jointly represented by fW) with the sleep- 
promoting ventrolateral preoptic nucleus (VLPO, /S) [17]. The 
model VLPO population expresses the inhibitory neurotransmitter 
GABA [cS) while the neurotransmitter expressed by the model 
wake population represents the joint effects of the transmitters 
expressed by LC, DR and TMN, namely norepinephrine, 
serotonin and histamine, respectively (cW). We constructed the 
model using our previously developed neuronal population firing 
rate and neurotransmitter formalism [24] (see Methods and Model 
section). Briefly, in this formahsm, the firing rate of a pre-synaptic 
population, /Jr(f) (in Hz), induces expression of neurotransmitter 
concentration, cX{f), which, in turn, acts as input to post-synaptic 
populations (see Eqs. (1), (2) with X,Y = W ,S). 

In the flip-flop, transitions between sleep and wake states are 
governed by the homeostatic sleep drive, h{t), that describes the 
universally recognized propensity for sleep that increases during 
time awake and decreases during sleep, and is thought to involve 
the neuromodulator adenosine (reviewed in [29]). As such, h 
increases when fW is at a high level [fW >6w) simulating the 
wake state (Fig. IB, see Eq. (5)). Increasing values of A cause the 
sleep-promoting population to activate and force a transition to 
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Figure 1. A: Schematic for the single wake/sleep flip-flop model. Rectangles represent neuronal groups and are labeled with their firing rate 
variables: fW for wake-promoting and fS for sleep-promoting. Circles represent neurotransmitter concentrations expressed by the neuronal groups 
and are labeled with the variable names: c W for the monoaminergic transmitters of the wake-promoting populations and cS for GABA expressed by 
the sleep-promoting population. The triangle represents the homeostatic sleep drive (/?). Solid lines ending in filled circles represent inhibitory 
synaptic projections. Large open arrow (S(t)) represents random excitatory stimuli to the wake-promoting population from sources external to the 
network. B: Time traces for fWJ'S, cW, cS (upper panel) and /; (lower panel). C: Bifurcation diagram of fixed point solutions of the fast subsystem 
(Eqs. (1), (2) with h treated as a fixed parameter. The upper (red) and lower (blue) branches indicate stable solutions which are joined through a 
branch of unstable solutions (gray dashed) at saddle-node bifurcations occurring at b = RKtY (right knee) and h = LKw (left knee). The orange curve 
is the trajectory of the full model projected onto the h—fW plane. Horizontal dashed line indicates Ow, the fW threshold defining the wake state. 
doi:1 0.1 371 /journal.pone.0094481 .gOOl 



the state in which fS is at a high level and fW is at a low level, 
simulating sleep. As fW drops to a low level ifW <dw), h 
decreases until it deactivates fS causing a transition back to the 
fW dominant state or simulated waking. 

Relative to the time scale of transitions between the simulated 
waking and sleep states, the homeostat h(t) is slowly varying. For 
this analysis, we consider a separation of timescales, in which the 
variables fW, fS, c W and cS vary on a timescale faster than that 
of A. We call Eqs. (1) and (2) governing variables (fW /S,cW,cS) 
the fast subsystem, and Eq. (5) governing h the slow subsystem. 
Here, we briefly summarize the intuition behind this fast-slow 
decomposition; for a rigorous review see [30-32] and for a 
previous analysis on a related model see [33]. 

We treat A as a slowly varying parameter of the fast subsystem. 
For a fixed value of h, solutions to the fast subsystem approach a 
stable fixed point, the value of which may depend on the initial 
condition. The value of fW at this fixed point dictates if h will 
increase or decrease in the fuU model. Slow variations in h are 
instantly reflected by the convergence of the fast subsystem to a 
new stable fixed point. Thus, we may track the trajectory of the fuU 
model as a slow evolution through stable fixed point solutions of 
the fast subsystem. 

The fixed point solutions, for fW, of the fast subsystem as a 
function of h form a Z-shaped curve, as plotted in the bifurcation 



diagram in Fig. 1 C . The upper (red) and lower (blue) branches are 
stable solutions. They are joined by an unstable branch of 
solutions (dashed) at saddle-node bifurcations occurring at a high 
value oi h, h = RK\Y referred to as the "right knee", and at a low 
value oih, h = LKfff referred to as the "left knee". Plotted on top 
of these fixed point solutions of the fast subsystem is the trajectory 
of the fuU model when h is allowed to evolve according to the slow 
subsystem (light blue curve). The Z-shaped curve of fast subsystem 
solutions forms the basis for hysteresis loop dynamics of the fuU 
model. When fW is at a high level (ffV >Bw, dotted line), h slowly 
increases and the fast sub.system variables remain attracted to the 
upper branch of fixed points. As h increases beyond the right knee 
of the curve at h = RKw, the only stable solution of the fast 
subsystem corresponds to low values of fW and the trajectory 
quickly approaches the lower branch of fixed points. On this 
branch fW<6w, so h decreases and the full model trajectory 
tracks the lower branch of frxed points until the saddle-node 
bifurcation at h = LK\Y- As h decreases below the left knee, the 
only stable solutions for the fast subsystem are on the upper branch 
and the trajectory jumps up, thus completing one cycle of the 
hysteresis loop, or one sleep-wake cycle. 

Bout durations and the influence of variability. In the 
flip-flop model during one sleep-wake cycle, or one cycle around 
the hysteresis loop, wake and sleep bout durations are determined 
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Figure 2. Effects of model variability on hysteresis loop dynamics of a single flip-flop. A, B: Variable neurotransmitter expression levels 
modulate position of the knees of the Z-shaped curve of fixed point solutions of the deterministic model (shown for reference with the deterministic 
trajectory (orange curve) in A). The projection of the noisy model trajectory (A: black curve) onto the bifurcation diagram illustrates the perturbations 
to the underlying hysteresis loop that result in variable wake and sleep bout durations illustrated in the time traces for fW,fS, cW, cS (B: upper 
panel) and /; (B: lower panel). C,D: Random excitatory inputs to the wake-promoting population perturb the model trajectory (C: black curve) as it 
evolves around the deterministic hysteresis loop (Z-shaped curve of fixed point solutions and deterministic trajectory (orange curve) shown for 
reference in C), resulting in wake bouts with brief and longer durations, and sleep bouts with variable durations, as illustrated in the time traces for 
fW,fS, cW, cS (D: upper panel) and /; (D: lower panel). 
doi:10.1371/journal.pone.0094481.g002 



by the time the trajectory takes to traverse the upper and lower 
branches, respectively, of the Z-shaped curve of fixed point 
solutions of the fast subsystem (Fig. IC). Thus, bout durations are 
governed by the distance \RKw — LKw\ and the time dynamics of 
h. Asymmetry in wake and sleep bout durations can be introduced 
by different values of h time constants during its increasing and 
decreasing phases, dictated by T/,_„^ and Xh,doKn in our model. The 
time evolution of A, and thus bout durations, are further influenced 
by the values of RKw and LKw relative to the maximum and 
minimum limits on values of /7, huf^x and hmin, respectively. Since h 
follows exponential dynamics, the evolution of h is slower when h 
approaches either h^ax or h,„i„, and is faster if h is near h„nn but 
increasing or h is near h^ax but decreasing. Thus, asymmetry in 
bout durations is also introduced due to the location of the Z- 
shaped curve within the interval {hminMmax)- 

As an example of how these two factors can compete in 
introducing asymmetry of bout durations, the longer wake bouts 
(460s) compared to sleep bouts (280s) shown in Fig. IB were 
obtained with a faster h time constant during simulated wake 
{T^h.up = 600s) than during simulated sleep {zh,down = 700s), but with 
a shorter distance between RKw and h„,ax than the distance 



between LKw and . Hence, despite a faster time constant for h 
while the trajectory is on the upper branch of the Z-shaped fixed 
point curve, the evolution of h is slower than on the lower branch 
because of the proximity of A values to h„,ax- 

A physiological flip-flop switch would be subject to different 
sources of variability. Here, we summarize our analysis of how 
physiologically motivated sources of noise, which we include in our 
models [24], perturb hysteresis loop dynamics; see Appendix SI 
for complete details. First, we consider the effects of variability in 
neurotransmitter expression, modeled by multiplicatively scaling 
the steady state neurotransmitter expression functions cX^ 
[X = W ,S) in Eq. (4) by the randomly varying term ^x(0- The 
amplitude of Itx randomly varied (with uniform distribution and 
unit mean) at discrete times dictated by a Poisson process. In the 
single wake/sleep flip-flop model, (£,s) affects the steady state 
expression levels of monoamines (GABA) by the wake-promoting 
(sleep-promoting) population and thus modulates the level of 
synaptic inhibition between populations. As £,w, is take on 
different values, in an interval around 1, they affect the Z-shaped 
frxed point curve, and thus the hysteresis loop, by changing the 
values oi RKw and LKw- The distance between RKw and LKw 
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Figure 3. Hypnograms (A,B) and summary statistics (C) of experimental rat sleep recordings and model simulations in baseline 
sleep conditions and after 24 h of REM sleep deprivation. A,B: Representative experimental and simulated hypnograms depicting transitions 
between states of wake, NREM sleep and REM sleep in the baseline condition (A) and after REM sleep deprivation (B). C: Summary statistics, such as 
mean bout durations, mean number of bouts and mean percent time spent in state, for experimental data (unhatched bars) and model simulations 
(10 runs for 4 simulated h's, hatched bars) in the baseline (white bars) and post REM sleep deprivation (gray bars) conditions. No significant 
differences were computed between data and model results (paired, 2-taiied t-tests, /)<0.05). 
doi:10.1371/journal.pone.0094481.g003 



decreases for and values less than 1 and increases for values 
greater than 1 (Fig. SI). This effect on the width of the hysteresis 
loop i.s a direct result of decreases and increases in inhibition 
between populations. 

As (^f^ and ^5 independently vary randomly around 1, the 
hysteresis loop randomly changes shape and position leading to 
variations in wake and sleep bout durations (Fig. 2A,B). Given that 
variations m^^r and induce both lengthening and shortening of 
the hysteresis loop, one might expect that the variable bout 
durations would be symm(;tri( ally distributed about the bout 
durations dictated by the deterministic model. Our analysis 
indicates, however, that the distributions of bout durations depend 
sensitively on the relative time scales of the variability (i.e. the 
average frequency oi £,w and random variations) and the 
deterministic bout durations. In the noisy flip-flop model, the 
majority of bout durations are shorter than the durations of the 
deterministic model and they are distributed with a tail of longer 



durations (Fig. S2). Differences in the location of the Z-shaped 
fixed point curve in the {h„i„,h„ax) interval can introduce 
significant differences in the extent of the positive tail for wake 
or sleep bouts. 

The second source of physiologically motivated variability we 
include in our model is brief, random excitatory stimuli input to 
the wake-promoting population, 5(1), representing external inputs 
from brain areas not included in the model such as sensory or 
cortical areas. The inputs have random amplitude, decay with a 
fixed time constant and occur at discrete times dictated by a 
Poisson process. They do not modulate the hysteresis loop, but 
instead perturb the trajectory as it evolves around the deterministic 
hysteresis loop (Fig. 2G,D). If a stimulus occurs during a wake 
bout, the trajectory is perturbed to higher fW values which may 
extend the wake bout if the trajectory is close to RKw- If the 
stimulus occurs during a sleep bout, it may result in a brief 
activation of fW with a return to the sleep state, which we define 
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Table 1. Probabilities of behavioral state transitions. 






Data 




Model 






Baseline 


Post REM-dep 


Baseline 


Post REM-dep 


From Wake To NREM 


0.9742 ±0.0577 


0.7358±0.1405 


0.9816±0.0311 


0.7950 ±0.0877 


To REM 


0.0258 ±0.0577 


0.2642 ±0.1405 


0.0184±0.0311 


0.2050 ±0.0877 


From NREM To Wake 


0.6449±0.0821 


0.2331 ±0.0561 


0.6705 ±0.0985 


0.2154±0.0679 


To REM 


0.3551 ±0.0821 


0.7669±0.0561 


0.3295 ±0.0985 


0.7846 ±0.0679 


From REM To Wake 


0.5962±0.1927 


0.6236±0.1894 


0.6246 ±0.1543 


0.5190±0.0802 


To NREM 


0.4038±0.1927 


0.3764±0.1894 


0.3754±0.1543 


0.48 10 ±0.0802 


Probabilities of behavioral state transitions computed from the experimental sleep recordings under baseline and post-REM 
computed from 10 simulation runs of the coupled flip-flop model under each simulated condition. 


sleep deprivation conditions, and 
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as a brief wake bout, or result in a full transition to the wake state. 
Which result occurs is influenced by the position of the trajectory 
on the lower branch of the Z-shaped fixed point curve when the 
stimulus arrives. StimuK occurring when the trajectory is closer to 
LKw usually result in a transition to wake and stimuli occurring 
when the trajectory is closer to RKw usually result in a brief wake 
bout. This dependence is due to whether the stimulus pushes fW 
above the middle branch of unstable fixed point solutions which 
generally acts as a separatrix between the basins of attraction of 
the upper and lower branches of stable frxed point solutions. 

The random excitatory stimuli to the wake-promoting popula- 
tion have different effects on the wake and sleep bout duration 
distributions (Fig. S2). The wake bout distribution is bimodal with 
a peak at very short durations indicating the brief wake bouts 
induced by the stimuli and a peak at longer durations for the bouts 
generated through hysteresis loop dynamics. The distribution of 
sleep bout durations takes on a more exponential-shape because 
sleep bouts are interrupted at random times by either brief wake 
bouts or early transitions to wake. Thus, an external input that 
preferentially targets one of the flip-flop populations can yield 
sleep-wake patterns with qualitative differences in wake and sleep 
bout durations. 

Coupled flip-flop model for rat sleep-wake regulation 

As the above analysis suggests, a single flip-flop with multiple 
sources of variability and appropriately tuned parameters would 
be able to generate transition dynamics between states of wake and 
sleep, or between states of NREM and REM sleep, similar to those 
observed in the rat. The goal of this study was to investigate how a 
wake/sleep flip-flop and a REM-on/REM-oflF flip-flop can be 
coupled to reproduce transition dynamics among the three states 
of wake, NREM sleep and REM sleep in rat. To constrain model 
parameter settings and network structure, we analyze experimen- 
tal recordings of rat sleep behavior under ad libitum (baseline) 
conditions and during REM sleep rebound after 24 h of REM 
sleep deprivation. Specifically, we consider data collected during a 
4 h window in the light period (rest phase) at the same circadian 
phase for both conditions. We assume minimal modulation of 
sleep-wake behavior by the 24 h circadian rhythm during this 4 h 
window and, thus, do not include the influence of the circadian 
rhythm in the model. As we describe below, analysis of the data 
motivated a minimal set of projections between the flip-flops to 
account for sleep-wake patterns in both conditions. These 
projections include an effect of activity of the wake population 
on the REM sleep homeostat and an excitatory effect of activation 
of the REM-on population to the wake-promoting population. 



As shown in Fig. 3, simulation results of our coupled flip-flop 
model were able to reproduce sleep-wake dynamics that were 
statistically similar to the experimentally recorded rat behavior 
under both conditions. Specifically, mean bout durations, number 
of bouts and percent time spent in each state for wake, NREM 
sleep and REM sleep, as well as probabilities (Table 1) were similar 
to the experimental data (paired, 2-tailed t-tests, ^<0.05). In this 
section, we first describe the key features of the experimental data 
that informed the model network structure and parameter settings. 
We then describe how the model was constructed to account for 
these key features, resulting in model dynamics similar to the data. 

Rat sleep-wake behavior in baseline conditions and after 
24 h REM sleep deprivation. As previously reported for these 
experiments [34], rats showed a statistically significant increase in 
the percent time spent in REM sleep after 24 h of REM sleep 
deprivation compared to baseline conditions (Fig. 3C). Additional 
analysis revealed that the REM sleep increase was due to a large 
increase in the number of REM bouts as well as an increase in 
mean REM bout durations. The percent time spent awake also 
significantiy decreased in the post-REM deprivation condition due 
to a trend towards shorter wake bouts. In particular, in baseline 
sleep conditions each rat had at least 1 wake bout of duration 
around 30 minutes, and 2 of the 5 rats had maximum wake bout 
durations of an hour or longer, accounting for the large variability 
in mean wake bout durations. In contrast, after REM sleep 
deprivation the longest wake bout any of the 5 rats exhibited was 
just over 30 minutes and the variability in durations was much 
reduced. This higher pressure for sleep in general and REM sleep 
particularly suggests that the REM sleep deprivation did not 
exclusively affect REM sleep homeostasis, but also influenced 
NREM sleep homeostasis. 

Analyzing the probabilities of behavioral state transitions from 
the experimental recordings under baseline and post-REM sleep 
deprivation identified some key features of state transition 
dynamics that we used to construct the interactions between the 
wake/sleep and REM-on/REM-oflF flip-flops in our model 
(Table 1). As is characteristic of normal rodent sleep patterning, 
the probability was higher for the termination of a REM sleep 
bout by a transition into the wake state than by a transition to 
NREM sleep. The data suggests that the responsible mechanisms 
must be fairly robust as the probabihties were very similar in the 
two conditions. In the baseline condition, the data exhibited the 
normal sleep pattern in which sleep initiates in NREM sleep with 
the transition to REM sleep occurring after a latency period. 
Specifically, 4 out of 5 rats always transitioned from wake to 
NREM sleep. For the 1 rat that showed several wake-REM 
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Figure 4. Schematic of the coupled flip-flop sleep-wake regulatory network model. Rectangles represent neuronal groups and are labeled 
with their firing rate variables: flV for wake-promoting, /S for sleep-promoting, /K"" for REM-promoting (REM-on), /K"'^ for REM-suppressing (REM- 
off). Circles represent neurotransmitter concentrations expressed by the neuronal groups and are labeled with the variable names: clV for the 
monoaminergic transmitters of the wake-promoting populations, cS, cR"" and cR"f^ for GABA expressed by the sleep-promoting, REM-promoting 
and REM-suppressing groups, respectively. Triangles represent the homeostatic drives for sleep (h) and REM sleep (stp). Solid lines ending in filled 
circles represent inhibitory synaptic projections; dashed lines indicate our proposed interactions between the flip-flops (filled circle indicates activity- 
suppressing, arrow indicates activity-promoting). Large open arrow (S(t)) represents random excitatory stimuli to the wake-promoting population 
from sources external to the network. B,C: Time traces for fW ,fS (B, upper panel), /; (B, lower panel), (C, upper panel) and stp (C, lower 
panel) for the simulated baseline sleep shown in Fig. 3A. 
doi:10.1371/journal.pone.0094481.g004 



transitions, these transitions occurred exclusively as an interrup- 
tion of REM sleep by a brief wake bout of duration 10-30 s, thus 
occurring as a REM-wake-REM transition. In the post-REM sleep 
deprivation condition, all 5 rats exhibited wake-REM transitions 
with the mean transition probability increasing to just over 25%. 
Again, these transitions occurred exclusively, in all rats, as REM- 
wake-REM transitions where the intervening wake bout was of 
brief duration (10-30 s). These findings suggest that the mecha- 
nisms by which REM sleep is terminated promote initiation of 
activity in wake-promoting populations but do not participate in 
the maintenance of that activity. 

Replicating stereotypical wake-NREM-REM transition 
dynamics. In this section we describe how obtaining the 
stereotypical state transition pattern of wake to NREM sleep to 
REM sleep motivated the inclusion of an effect of activity of the 



wake population on the REM sleep homeostat in our coupled flip- 
flop model. As a starting point, we describe the dynamics of the 
coupled model using fast-slow decomposition. The coupled flip- 
flop model (Fig. 4A) consists of two pairs of neuronal populations, 
each pair reciprocally coupled by neurotransmitter-mediated 
inhibition, representing a wake/sleep flip-flop, with population 
firing rates and neurotransmitter levels modeled by Eqs. (1) and (2) 
with X,Y=W,S, and a REM-on/REM-oflf flip-flop, with 
population firing rates and neurotransmitter levels modeled with 
Eqs. (1) and (2) with X ,Y = R"" ,R''ff . The REM-on population 
{fR°") represents the joint activity of identified REM sleep- 
promoting populations and its neurotransmitter [cR°") represents 
their joint GABAergic signaling. The REM-off population ifR"^^) 
and its neurotransmitter {cR°^^) represent the GABAergic signal- 
ling of the identified areas with REM-ofF activity. Transitions in 
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Figure 5. Coupled flip-flop model dynamics replicating wake-NREM-REM transition pattern. Hypnogram depicting state transition 
dynamics (A), homeostatic sleep drive /?(/) (B) and REM sleep homeostatic drive stp(t) (C) in the deterministic model (all sources of variability 
removed) during REM/NREM cycling (red portion of curves) followed by a prolonged wake bout (light blue portion) and a subsequent return to REM/ 
NREM cycling (orange and dark blue portion). D: Z-shaped curve oifW fixed point solutions (Z) of the fast subsystem of the wake/sleep flip-flop with 
h as a bifurcation parameter, with the projection of the trajectory of the full model onto the h-f]¥ plane. The dashed line corresponds to the 
threshold Om^ determining the wake state E: S-shaped curve of fixed point solutions (S) of the fast subsystem of the REM/NREM flip-flop, with stp as a 
bifurcation parameter, with the projection of the trajectory of the full model onto stp-fR"" space. 
doi:10.1371/journal.pone.0094481.g005 



the sleep-wake flip-flop are governed by the homeostatic sleep 
drive, h{t), modeled by Eq. (5). As described above, h increases 
during waking {fW>6w) to promote activation of the sleep- 
promoting population and the transition into sleep, and decreases 
during sleep ifW <6w) to promote deactivation of the sleep- 
promoting population. In the REM-on/REM-ofF flip-flop, tran- 
sitions are governed by the REM sleep homeostatic drive modeled 
by the variable stp, as a reference to the hypothesized short-term 
process involved in REM sleep homeostasis [15], using Eq. (6). To 
replicate the reported phenomena of REM sleep homeostasis, we 
model stp as increasing during NREM sleep (fR"" <Qr<>«) to 
promote deactivation of the REM-off population and the 
transition into REM sleep, and decreasing during REM sleep 
{/K"" > 6/;o«) to promote activation of the REM-off population. 
This implementation of the REM sleep homeostatic drive is 
consistent with the concept that NREM-REM cycling is a sleep- 
dependent process [9,15] and it generates cychng solutions that 
are robust to variations in parameter values, such as strength of 
inhibition between the REM-on and REM-off populations, as 
shown in previous analysis [27]. Additionally, we include 
variability in all neurotransmitter expression and brief random 
excitatory stimuli to the wake-promoting population ((5(?) in 
Fig. 4A). 



Qualitatively we can understand the dynamics of these two flip- 
flops during normal sleep behavior as follows (Fig. 4B,C): during 
the wake state, the wake-promoting population is activated and h 
increases. To mimic the activity of identified physiological REM- 
off neuronal areas [8] , the REM-off population is activated and 
the REM-on/REM-ofF flip-flop does not exhibit any transitions. 
When increasing h forces activation of the sleep-promoting 
population, a sleep episode is initiated, h starts to decrease and 
stp starts to increase. A sleep episode may be terminated at any 
time by deactivation of the sleep-promoting population as 
governed by h and a return to the wake state, or may be 
interrupted by a brief wake bout generated by the random 
excitatory stimuli to the wake-promoting population, &(t). As a 
sleep episode continues, however, increasing stp will deactivate the 
REM-off population, allowing the REM-on population to activate 
and a REM bout to occur. During a REM bout, stp decreases and 
the REM bout can terminate due to activation of the REM-off 
population which returns the model to the NREM sleep state. 
Alternatively, the REM bout can be terminated by a brief wake 
bout. Eventually, the sleep episode ends when h reaches a 
sufficiently low level resulting in deactivation of the sleep- 
promoting population and reactivation of the wake-promoting 
population. 
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To quantitatively understand transition dynamics in the coupled 
flip-flop model, we extend the fast-slow decomposition analysis to 
define two hysteresis loops that dictate trajectory dynamics, one 
defined by the wake/ sleep flip-flop and the other by the REM-on/ 
REM-ofT flip-flop. Informally, we consider both homeostatic drives 
as slow variables in the coupled flip-flop model with the 
homeostatic sleep drive h acting on a slower time scale than the 
REM sleep homeostatic drive stp. The firing rate and neuro- 
transmitter level variables for the wake- and sleep-promoting, 
REM-on and REM-off' populations {fW, cWJS, cS, fR°", cR"", 
compose the fast subsystem. Since wake/sleep and 
REM-on/REM-oflf transitions occur during distinct phases of the 
trajectory and since each homeostatic drive only affects the fast 
variables of one of the flip-flops, we can separately apply fast-slow 
decomposition to each flip-flop to define two hysteresis loops. For 
the wakc/sk'<;p flij)-fl()p, the fast-slow decomposition is identical to 
that described for the single flip-flop model, namely by considering 
h a frxed parameter in Eqs. (1) and (2) with X ,Y = IV, S, we obtain 
a Z-shaped curve of fW fixed point solutions. We refer to this 
curve as Z. For the REM-on/REM-oflf flip-flop, we consider stp a 
fixed parameter in Eqs. (1) and (2) with X,Y = R°",R''-0' and obtain 
an S-shaped curve of fixed point solutions of fR"", which we refer 
to as S. The reversal in the shape of the curve of fixed point 
solutions occurs because stp increases when fR"" is deactivated 
(during NREM sleep). Solution trajectories of the model can be 
tracked in relation to these two fixed point curves, Z and S, such 
that wake and sleep behavior corresponds to the trajectory' 
evolving along the upper and lower branches, respectively, of Z 
and, during sleep, REM and NREM episodes occur as the 
trajectory evolves along the upper and lower branches, respec- 
tively, of S. 

To achieve activation of the REM-off population during wake 
and to robustly generate the stereotypical sleep pattern of wake- 
NREM-REM transitions that occurs after an extended period of 
waking, we need to include in the model a mechanism that forces 
the trajectory to remain on the lower branch of 5 during wake and 
ensures that when the transition to sleep occurs, stp is at a value 
less than the right knee of 5. To determine these constraints on the 
model, we consider the physiological hypotheses for the action of 
the wake state on the REM sleep homeostatic drive. As discussed 
in [35,36], there is debate whether the REM homeostatic drive 
increases during wake, decreases during wake or is unaffected by 
the wake state. In the context of our coupled flip-flop model, we 
can implement these hypotheses through the effect on stp of 
activation of the wake population. If stp increases during wake, 
this may lead to an immediate (or almost immediate) transition 
from wake to REM, since stp may reach values close to or above 
the right knee of S during an extended wake bout. Similarly, stp 
remaining constant during wake could also lead to immediate (or 
almost immediate) transitions from wake to REM. For example, as 
the data show, the majority of REM bouts end in a transition to 
wake. Thus, at the transition to wake stp is at a level between the 
two knees of S, and could be at a value close to the right knee. 
Even if the trajectory is forced down to the lower branch of S 
during wake in order to activate the REM-ofF population, when 
the next sleep episode occurs there could be a very short transition 
to REM as stp could quickly evolve to the right knee of S. To 
guarantee a finite REM latency, as is stereotypical, we model that 
stp decreases during wake to a reset value stpr- This mechanism 
corresponds to a saturating inhibitory effect of the wake state on 
the REM sleep homeostatic drive and insures that the REM-off 
population is activated during wake. We model this mechanism as 
an additional condition on stp dynamics that when fW>9w 
(wake), stp decays to stpr with a time constant Xstpw> when 



fW <dw (sleep) stp dynamics is given by Eq. 6 (see Eq. S15 in 
Appendix SI). 

To illustrate how this condition on stp due to the wake state 
creates the stereotypical pattern of wake-NREM-REM transitions, 
we consider the trajectory of the deterministic model shown in 
Fig. 5 as a hypnogram (parameter values differ slightiy than those 
in Table SI). At the beginning of the simulation, the model is in 
the sleep state exhibiting regular transitions between NREM and 
REM sleep (A, red portion of the cur\'e). The trajectory^ is slowly 
evolving leftwards on the lower branch oi Z (D, red) as h decreases 
(B, red), and stp (C, red) is alternating between values associated 
with the left and right knees of S as the trajectory traverses the 
hysteresis loop defined by S (E, red). When h decreases below the 
left knee of Z, the wake/ sleep flip-flop tiansitions to the wake 
state, interrupting the current REM or NREM bout (all panels, 
light blue). The trajectory jumps to the upper branch of Z, h 
begins to increase and stp is driven down to stp^ {stpr = 0 in this 
simulation), forcing the trajectory on the lower branch of 5. At the 
end of the wake bout as h increases beyond the right knee otZ, the 
trajectory jumps down to the lower branch of Z, h starts of 
decrease, stp is released from its reset value stp, and begins to 
increase (all panels, orange). This portion of the trajectory 
represents the REM latency period. The first REM bout occurs 
when stp reaches the right knee of S and the trajectory begins 
traversing the hysteresis loop defined by S (all panels, dark blue). 
For values of stpr less than the left knee of S, the REM latency 
period can be of longer duration than the interval between NREM 
and REM bouts determined by the hysteresis loop dynamics of S. 

In the full model with variability sources included, this 
inhibitory effect of the wake state on the REM sleep homeostatic 
drive remains a robust mechanism promoting the wake-NREM- 
REM sleep transition pattern. In particular, we set stpr sufficientiy 
low [stpf = 0) so that it remains less than the right knee of S despite 
the modulation of the knees of S induced by the neurotransmitter 
variability. 

There are other possible ways to robustiy generate the Wake- 
NREM-REM transitions achieved here, such as including 
projections from the sleep-promoting population to the REM-on 
or REM-off populations with a decaying inhibitory effect or a 
decaying excitatory effect, respectively, that would prevent 
activation of the REM-on population at the onset of a sleep 
episode regardless of the influence ot stp. However, the inclusion 
of such additional projections would increase the complexity of the 
model and introduce new timescales and other parameters. In the 
absence of experimental hypotheses supporting such additional 
mechanisms, we focused on the effect of the waking state on the 
REM sleep homeostatic drive to generate the appropriate state 
transition dynamics. 

Generating REM-wake transitions. Behavioral state tran- 
sition probabilities computed from the experimental sleep 
recordings suggest a robust mechanism that terminates REM 
bouts with a transition to waking instead of NREM sleep. In the 
model, two mechanisms can generate transitions to the wake state 
from either the NREM or REM sleep state: either h decays below 
the left knee of Z forcing the activation of fW or an input to the 
wake-promoting population perturbs fW off the lower branch of 
Z causing its activation and disruption of the sleep state. If these 
mechanisms are independent of the activity of the REM-on and 
REM-off populations, we expect that most sleep-wake transitions 
would occur from the NREM state because the model spends 
much more time in NREM compared to REM sleep. To generate 
REM-wake transitions, then, there could be an excitatory effect of 
activation of the REM-on population to the wake-promoting 
population, or an inhibitory effect to the sleep-promoting 
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population. Within the context of our coupled flip-flop model, 
these excitatory and inhibitory effects could be direct synaptic 
projections from the REM-on population to the wake- or sleep- 
promoting population, REM-dependent actions on the homeo- 
static sleep drive or REM-dependent actions that indirectiy 
influence transient or sustained external inputs to the wake- or 
sleep-promoting populations. 

Investigation of these different effects in the model indicated 
that direct synaptic projections from the REM-on population to 
either the wake- or sleep-promoting population, or actions on the 
homeostatic sleep drive promoted the occurrence of wake 
transitions by modulating the hysteresis loop Z or h dynamics 
such that h decreased below the left knee oi Z, thereby causing the 
trajectory to jump to the upper branch of Z. Such transitions from 
sleep to wake almost always led to an extended wake bout. 
Including a REM-dependent external input to either the wake- or 
sleep-promoting populations whose activity was sustained during 
REM bouts similarly affected model dynamics to cause transitions 
to extended wake bouts. While these mechanisms did not affect the 
ability of the model to replicate the experimental data in the 
baseline condition, it prevented replication of the post-REM sleep 
deprivation condition where REM-wake-REM transitions with 
brief intervening wake bouts were prevalent. 

To investigate the effects of transient REM-dependent external 
inputs, we modulated the random excitatory stimuli to the wake- 
promoting population, 5{{), dependent on activation of the REM- 
on population. We found that model solutions fit the statistics of 
the experimental sleep-wake patterns in both baseline and post- 
REM sleep deprivation conditions when the frequency of &(t) 
stimuli increased during the REM sleep state. Hence, the high 
REM-wake transition probability and the occurrence of REM- 
wake-REM transitions in the data indicated an indirect effect of 
REM-on activity on the wake /sleep flip-flop that initiated 
activation of the wake-promoting population but not maintenance 
of its activity. To implement this mechanism, we could include an 
additional source of transient excitatory stimuli to the wake 
promoting population that is only active during REM sleep. 
However, in the interest of keeping the model compact, we 
modeled the interaction through our existing mechanism 5(t). See 
the Discussion for possible physiological motivations for such 
stimuh. Alternatively, including a source of REM-dependent 
transient inhibitory inputs to the sleep-promoting population 
similarly resulted in high probabilities of REM-wake transitions 
and the occurrence of REM-wake-REM transitions but the 
mechanism was less robust as it depended on sufficient inactivation 
of the sleep-promoting population that could result in sufiicient 
activation of the wake-promoting population. 

We modeled the effect of REM-on activation on the random 
excitatory stimuli to the wake-promoting population, d{{), with an 
increase in the frequency of occurrence of the stimuli during REM 
sleep. In this way, although the model spends less time in REM 
sleep, more frequent stimuli to fW during REM induce more 
REM-wake transitions without introducing excessive brief wake 
bouts during NREM. When an excitatory stimulus to the wake- 
promoting population occurs during a REM bout, fW is 
perturbed off the lower branch of Z. If h is sufficiently low, fW 
can transition to the upper branch of Z, which forces stp down to 
stp, terminating REM-on activity and the wake bout will be 
maintained as h evolves along the upper branch of .E. If h is closer 
to the right knee of Z when the stimulus arrives during the REM 
hout, fW may be briefly perturbed off the lower branch of .E but it 
will return to the lower branch, thus generating only a brief wake 
bout. The sleep state that the model returns to depends on the rate 
of decrease oistp to stpr induced by the brief /TF activation. If the 



rate of stp reset is fast, the model would return to the NREM sleep 
state because the trajectory would be forced below the left knee of 
S and fR"" would deactivate during the brief wake. If the rate of 
stp reset is slower, a REM-wake-REM transition can occur since 
stp would have remained above the left knee of S with perhaps 
only a transient reduction in REM-on activity. 

Accoimting for high variability of wake bout 
durations. In our experimental sleep recordings, especially 
under baseline conditions, some animals exhibited very long wake 
bouts. As described for the single flip-flop model, long wake bout 
durations can be obtained by adjusting time constants of h or 
providing the appropriate relationship between the saturation level 
of h and the knees of the Z. With variability of neurotransmitter 
expression included in the coupled flip-flop model, we can 
introduce very long wake bouts by setting the maximum saturation 
level of/), h„,i,\, below the h value of the right knee of Z. Without 
the neurotransmitter expression variability, these parameter 
settings would result in the model getting stuck in the wake state, 
because the model solution with h = hmax would be a stable fixed 
point of the fiall model. The neurotransmitter expression variability 
modulates the h value of the right knee oiZ and can move it below 
hmax to induce a transition to the sleep state. Thus, in this 
parameter regime, wake bout durations are overall longer than 
sleep bouts, they have higher variability and very long wake bouts 
are possible. In the Discussion we provide possible physiological 
mechanisms that could support these model dynamics. 

Differences between baseline and post-REM sleep 
deprivation sleep patterning. As shown in Fig. 3, the coupled 
flip-flop model is able to replicate sleep-wake patterning under 
both baseline and post-REM sleep deprivation conditions. The 
significant difierences in patterning between the two conditions 
include increases in percent time spent in REM sleep, the number 
and duration of REM sleep bouts and a decrease in the percent 
time spent in waking (2-tailed paired t-test, p < 0.05). We captured 
these differences in the model by adjusting parameters governing 
three model components. First, we increased the number and 
duration of REM bouts by adjusting the maximum and minimum 
saturation levels of stp. The maximum saturation level, stpmax> was 
increased to promote shorter latencies to REM-on activation and 
the minimum level, stp„,j„, was increased to slow down the 
trajectory's evolution on the upper branch of S during a REM 
bout. To further aflFect REM bout durations, we lengthened the 
REM-NREM hysteresis loop by increasing the distance between 
the knees of S. In our model formalism, the distance between the 
knees of S can be modified in different ways, including the 
addition of external inputs to either the REM-on or REM-off 
populations. Alternatively, we adjusted the influence of stp on 
/R"-^ activity (see Appendix SI for details). An additional 
consequence of a longer REM-NREM hysteresis loop is an 
increase in REM-wake-REM transitions, which also was a feature 
of post-REM sleep deprivation sleep patterning. As described 
above, REM-wake-REM transitions occur in the model when a 
random excitatory stimulus arrives to the wake-promoting 
population during a REM bout. As a result of brief fW activation, 
stp is forced to decrease towards the reset value stpr. However, 
when the knees of S are further apart, the decrease in stp is less 
likely to push stp below the left knee of S during the brief wake 
bout. When the brief wake bout ends, the REM-on/REM-ofi" flip- 
flop remains in the REM state and REM bouts are more robust to 
this kind of interruption. 

Secondly, we obtained the decrease in percent time spent in 
waking as a result of REM sleep deprivation through a decrease in 
wake bout durations. Specifically, we increased the maximum 
saturation level ofh, h„utx> to a value above the right knee of Z. As 
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described above, this shortened wake bouts by aUowing h to freely 
evolve past the right knee of Z to induce the deactivation of fW 
and end the wake bout. The third model component we adjusted 
to capture post-REM sleep deprivation patterning was to reduce 
the frequency of the random excitatory stimuli to the wake- 
promoting population. A large reduction in stimulus frequency 
during the wake and NREM states resulted in far fewer brief wake 
bouts during NREM sleep that fragmented NREM and prevented 
the increase of the REM sleep homeostatic drive stp. Therefore, 
we obtained the large increase in the NREM to REM transition 
probability exhibited in the sleep recordings (Table 1). The 
reduction in stimulus frequency during REM sleep, while not as 
great, reduced the interruption of REM bouts by terminating 
stimuh. In the Discussion, we provide possible neural mechanisms 
that could account for these parameter changes. 

Discussion 

The proposed conceptual models for the control of REM sleep 
by mutually inhibitory networks of REM-on and REM-olf 
populations leave a number of questions unanswered, particularly 
regarding the interactions between the populations controUing 
REM sleep and those controlling wake and NREM sleep [5,8,19]. 
In an attempt to shed some light on these questions, we 
constructed a simplified, yet cohesive, mathematical model based 
on mutually inhibitory flip-flop networks for the control of state 
transitions between wake, NREM and REM sleep states. We 
utilized a modeling formalism that correlates with the physiolog- 
ical network structure of neurotransmitter-mediated interactions 
among state-promoting neuronal populations. The network 
structure, in particular the interactions between the wake/sleep 
flip-flop and the REM-on/REM-ofl" flip-flop, was motivated by the 
behavioral state transitions observed in experimental rat sleep 
recordings. A minimal set of interactions was identified by 
constraining the network model to accurately replicate experi- 
mental state transition dynamics in both baseline sleep and post- 
REM sleep deprivation recovery sleep in a consistent manner. As 
discussed below, these modeled interactions predict physiological 
mechanisms that can be targeted in future experimental studies to 
more definitively address some of the unanswered questions of 
REM sleep regulation. 

A strength of this study was the use of experimental rat sleep 
recordings to motivate the construction of the network model. We 
relied on replication of summary bout statistics and probabilities of 
behavioral state transitions to constrain the model. Recently, 
higher order statistical approaches, such as survival-based analysis 
of wake and sleep bout durations, have been applied to sleep 
recordings to identify effects of disease states and experimental 
manipulations [37-40]. To compare survival analyses of bout 
durations of the data and model results requires a sufficient 
number of bouts [41]. Since our data sets contained only five 4 h 
recordings in each condition, baseline and post-REM sleep 
deprivation, there were insufficient numbers of bouts, particularly 
REM sleep bouts, for a bout duration survival analysis to be 
meaningful. We note, though, that the standard summary statistics 
and state transition probabilities were sufficient to rule out possible 
interactions between the wake/ sleep and REM-on/REM-off flip- 
flops. Specifically, to robustly generate transitions from REM sleep 
to wake, an alternative interaction between flip-flops is a direct 
excitatory projection from the REM-on population to the wake 
population. However, simulations with this alternative model 
structure were not able to replicate all summary statistics and 
behavioral state transition probabilities. In particular, this alter- 
native model could not simultaneously replicate the correct REM- 



wake transition probability and number of REM bouts in either 
the baseline or post-REM sleep deprivation condition. In this 
alternative model, REM-wake transitions most often resulted in an 
extended wake bout instead of a brief wake bout as part of a REM- 
wake-REM transition that occurred more often in the data. Thus, 
we believe that summary statistics and transition probabilities were 
adequate to constrain our simplified model. 

Mathematically, a mutually inhibitory flip-flop network pos- 
sesses the dynamics of a hysteresis loop. The inherent symmetry of 
a hysteresis loop and the regularit}' of its dynamics may call into 
question its suitability for replicating the highly variable state 
transitions of rat sleep-wake beha\ior. For example, survival 
analysis of wake and sleep (NREM and REM sleep combined) 
bout durations of rodent sleep revealed an asymmetry between 
wake and sleep bouts such that wake bout durations displayed 
power-law like distributions while sleep bout durations exhibited 
exponential distributions [28]. Such a qualitative difference in 
bout duration distributions suggests that wake and sleep state 
transitions are influenced by different mechanisms. Our analysis of 
the dynamics of a single ffip-flop network with physiologically 
motivated sources of variability included suggests that it can 
generate sufficient variability and asymmetry in wake and sleep 
bout duration distributions. In particular, an extended tail in wake 
bout distributions reffecting low numbers of very long wake bouts 
could be introduced by the appropriate relationship between the 
underlying hysteresis loop and the saturation level of the 
homeostatic sleep drive. Additionally, brief excitatory' inputs to 
the wake population generated a bimodal wake bout distribution 
and an exponential-like sleep bout distribution. Recent work has 
suggested that wake distributions may follow a multi-exponential 
distribution rather than a strict power-law distribution [42]. One 
can imagine that in the proper parameter regime, a bimodal wake 
bout distribution with an appropriate tail, as we have shown a 
single ffip-ffop with noise sources can generate, may result in a 
multi-exponential-like distribution. Further work fitting flip-flop 
models to survival analysis of experimental sleep data is clearly 
needed to assess the capabiHty of this network structure to account 
for all aspects of sleep state transition dynamics. 

Comparison to previous flip-flop models of the sleep- 
wake regulatory network 

Previous modeling studies have investigated a coupled flip-flop 
network for sleep-wake regulation in humans [23,25]. Our study 
differs in that we focus on accounting for the variable and 
polyphasic sleep-wake activity typical of rats rather than a 
caricature of stereotypical human sleep-wake activity. A more 
important difference involves the modeling formalism. The 
previous studies modeled the mean activity or firing rate of each 
neuronal population with the Morris-Lecar model, a generic 
model for single cell neuronal membrane potential [31,43]. When 
parameters are tuned appropriately, the Morris-Lecar model can 
generate fast transitions between states of low and high activity 
that is appropriate for modeling the activation and deactivation 
transitions of mean activity in neuronal populations. As a model 
for neuronal action potential generation, the Morris-Lecar model 
has the additional capability of generating osciUatory solutions. In 
the translation of using the model for population activity, this 
means that the model can generate spontaneous, regular 
transitions between low and high activity levels, without any 
external stimulation. Both the Rempe et al [23] and Kumar et al 
[25] studies exploit the model's capability of intrinsic population 
oscillatory dynamics to account for NREM-REM transitions. In 
the modeling formalism that we use, on the other hand, individual 
cell groups do not have the capability of intrinsic oscillatory 
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behavior. Transitions between low and high activity levels depend 
on changes in external inputs, specifically homeostatic drives for 
NREM or REM sleep. Thus, in our model results, NREM-REM 
transitions occur in response to changes in explicitly defined model 
components, instead of assumed population properties. As the aim 
of our study is to propose potential physiological interactions 
among behavioral state-promoting neuronal groups, we beheve 
that it is important to avoid incorporating extraneous assumptions 
into our modeling formalism. 

In these previous models, the coupling projections between the 
wake/ sleep and REM-on/REM-off flip-flops were similarly 
motivated by generating appropriate state transition dynamics, 
particularly the stereotypical transition pattern of wake to NREM 
sleep to REM sleep. The Kumar et al [25] model primarily 
achieves this pattern by the action of a REM sleep homeostatic 
drive whose dynamics were sensitive to wake and sleep states. This 
drive increased excitation to the REM-on population during sleep 
to promote its activation and decreased excitation during wake, 
which allowed the REM-off population to dominate at wake to 
sleep transitions. The Rempe et al [23] model included an indirect 
projection from the wake/sleep to the REM-on/REM-off flip-flop 
through an intermediary population, the extended VLPO 
(eVLPO). Wake and sleep dependent activity of the eVLPO gated 
activity of the REM-off population during wake and sleep states in 
order to promote its initial activation at the wake to sleep 
transition. Both models include additional feedforward connec- 
tions between flip-flops, such as direct inhibitory projections from 
the wake-promoting population to the REM-on population, that 
are suggested by anatomical studies and work to suppress REM-on 
activity during wake. We note that inclusion of similar additional 
feedforward projections from the wake/ sleep flip-flop to the REM- 
on/REM-ofF flip-flop in our model would have similar effects and 
not qualitatively change our results. In both the Rempe et al and 
Kumar et al models, replicating human sleep-wake patterns did 
not constrain the form of feedback projections from the REM-on/ 
REM-off flip-flop to the wake/sleep flip-flop. While an indirect 
feedback projection originating in the REM-ofiF population was 
included in the Kumar et al model, it was not necessary to obtain 
appropriate model behavior, and the Rempe et al model did not 
include any feedback projections. 

Since the Rempe et al [23] and Kumar et al [25] models 
simulated human sleep, they both included input from a circadian 
oscillator to contribute to the 24 h modulation of sleep-wake 
behavior. In this study, we focused on replicating variable and 
polyphasic rat sleep-wake behavior during a 4 h window in the 
rest phase, assuming minimal modulation by the circadian 
rhythm. Our modeling formalism, however, can readily include 
a neuronal population representing the suprachiasmatic nucleus 
(SCN) whose activity is driven by a circadian oscillator and which 
is coupled to the sleep- and wake-promoting populations within 
the network. Indeed, our previous work modeling rat sleep-wake 
behavior, in which REM sleep is generated by a reciprocal 
interaction network, suggested that multiple signaling pathways 
between the SCN and sleep-wake centers may be necessary to 
account for circadian modulation of rat sleep [44] . As a direction 
for future work, incorporation of the SCN circadian signal in the 
coupled flip-flop model to account for the 24 h variation of rat 
sleep-wake behavior would test the minimal interactions between 
the wake/ sleep and REM-on/REM-off flip-flops proposed here 
and perhaps identify additional constraints on the network 
structure. 



Model predictions 

To obtain robust wake-NREM-REM transition dynamics, our 
simplified coupled flip-flop network model predicts that during 
waking, the REM sleep homeostatic drive resets to a level 
corresponding to low REM sleep pressure. This prediction is 
similar to that of Benington and Heller [35] that REM sleep need 
is homeostatically related to NREM sleep rather than waking, 
such that in normal conditions it accrues during NREM and not 
waking. While Benington and Heller do not offer a mechanism 
defining the behavior of the REM sleep homeostat during the 
transition from sleep to wake, we found that providing a decay of 
the REM sleep homeostat to the reset level during wake best 
replicated rat sleep patterning. In the context of the coupled flip- 
flop model, this mechanism most parsimoniously pr()\'id<'s the 
appropriate dynamics to {'nsuri; that aftc'r an extended wake bout, 
the NREM sleep state is entered first before REM skx'p occurs. A 
physiological correlate of such a mechanism could be provided by 
the expression during sleep states of a substance mediating REM 
sleep need and its cessation during waking. The absence of 
expression during wake would lead to the degradation or uptake of 
the substance, decreasing its presence to low levels. Recent 
experimental results indicate that expression of MCH and 
Nesfatin-1, which are co-expressed l)y neurons in the tuberal 
hypothalamic area, exhibit sleep-dependent increases and wake- 
dependent decreases and have REM sleep-inducing effects 
[12,13,4.5]. Interestingly, while these substances are co-expressed, 
pharmacological experiments indicat(; that MCH promotes REM 
sleep [45] while Nesfatin-1 suppresses REM sleep [13]; however, 
recent optogenetic experiments report that acute activation of 
MCH neurons promotes maintenance of REM sleep [46]. An 
alternative physiological mechanism could be provided by sleep- 
dependent increases in activity of neurons that promote a 
transition into REM sleep. Neural activity in the median preoptic 
nucleus (MnPN) and the ventral lateral preoptic area (vlPOA) 
strongly correlates with REM sleep pressure [14]. These- lu-ural 
groups show higher activity during sleep than during waking 
which could mediate a wake-dependent decrease in REM sleep 
homeostatic drive. It is most likely that the physiological 
mechanism for a REM sleep homeostatic drive is more 
complicated than the simple, single drive variable included in 
the model. For example, experiments in cats have suggested that 
REM sleep pressure during REM deprivation and REM sleep 
rebound during the recovery from REM deprivation may be 
governed by mechanisms in separate brain areas [47]. Further 
experimental investigation can provide insight for a more accurate 
model. 

The robust propensity for REM bouts to terminate in a 
transition to wake exhibited in the data suggests the existence of an 
effect on wake-promoting populations by the activation of REM- 
on populations. However, while the anatomy has not been 
completely determined, there is little evidence of direct synaptic 
projections from the key identified REM-on populations to wake- 
promoting areas. For example, the efferent pathways of the SLD 
descend to areas that govern motoneuron activity to control 
muscle atonia during REM, and ascend to thalamic areas that 
induce REM cortical activation (reviewed in [8]). Thus, an indirect 
mechanism whereby wake-promoting populations are activated to 
terminate REM may be likely. We arrived at this conclusion by 
constraining the model to replicate transition dynamics of the 
data. As described above, model dynamics accurately replicated all 
summary statistics and probabilities of behavioral state transitions 
when activation of the REM-on population increased the external 
excitatory stimuli to the wake-promoting population, but not when 
it had a direct excitatory effect on the wake-promoting population. 
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A physiological correlate for the model mechanism could be 
provided by top-down projections to wake-promoting populations 
from cortical regions that are activated during REM sleep, but not 
NREM sleep. An alternative interpretation is that the wake- 
promoting population is more sensitive to external inputs during 
the REM sleep state, leading to a higher probability that wake- 
promoting populations will be briefly activated during REM sleep. 
Dynamically, increasing the rate of external excitatory stimuli 
during REM sleep is essentially equivalent. 

In the model, the occurrence of very long wake bouts and high 
variability of wake bout durations were achieved by allowing the 
homeostatic sleep drive h to reach its saturating limit hmax and 
relying on neurotransmitter expression variability to modulate the 
hysteresis loop to induce the transition out of the wake state. While 
this may not be a physiologically robust mechanism, it is, however, 
dynamically similar to the presence of a wake-promoting factor, 
such as orexin (see [48] for a review), that could disrupt or delay 
homeostatii:aily governed transitions to sleep. The time-depen- 
dence or variability in the expression of such a wake-promoting 
factor would then determine the transition out of sleep and thus 
the duration of the wake bout in a manner similar to random 
variations of neurotransmitter levels included in the model. 

The effects of 24 h REM sleep deprivation on rat sleep behavior 
were replicated in the model by modulating both the sleep and 
REM sleep homeostatic drives. Dynamics of the sleep homeostatic 
drive were modulated to promote transitions into sleep from 
waking. The need for this modulation may reflect a direct effect on 
sleep homeostasis by the REM deprivation protocol due to 
disturbances or slight losses of total sleep. For the effect on the 
REM sleep homeostatic drive, REM bouts were length(;n(;d and, 
consequendy, resistance to REM interruption was strengthened by 
modulating the influence of the REM homeostat on activation of 
the REM-off population. Such modulation could b(^ a mechanism 
for the proposed long-term process that rcgulat(;s tlu; daily amount 
of REM sleep [15]. Physiologically, it might reflect modulation of 
receptors for substances mediating short term REM homeostasis, 
such as MCH and Nesfatin-1, as a result of REM deprivation. 
This type of modulation may also be a mechanism to account for 
the resistance to REM interruption achieved by the reduction in 
the frequency of excitatory' stimuli to the wake-promoting 
population that was implemented to account for the effects of 
REM sleep deprivation. 

Conclusions 

In this study, we proposed a minimal model of an inhibition 
based network for the regulation of transition dynamics between 
the states of wake, NREM and REM sleep. We readily concede 
that the model may be too simple. Our intent, however, was to 
provide a cohesive inhibitory network structure, based on known 
physiology, and constrain the structure to account for experimen- 
tally measured sleep patterning in baseline conditions and after an 
experimental challenge, namely REM sleep deprivation. Key 
model results are predictions of the interactions between the 
subnetworks controlling sleep-wake transitions and REM-NREM 
transitions: feedforward effects from the wake/sleep subnetwork to 
the REM/NREM subnetwork that act to suppress REM 
propensity during waking, and feedback effects that indirectiy 
promote the initiation, but not maintenance, of waking as a result 
of REM sleep. The aim of this modeling study, similar to other 
recent physiologically-based modeling work on sleep-wake regu- 
latory networks [20-25], is to participate in the investigation of 
neuronal regulation of sleep and play the same role of providing a 
framework for understanding and interpreting experimental 
observations as the phenomenological two-process model [49] 



and the reciprocal interaction model [3] have done for the past 40 
years, but in the context of our increased knowledge of the 
underlying physiology of sleep-wake regulation. 

Methods and Model 

Experimental sleep recordings and sleep scoring 

Ethics statement. Experimental procedures were approved 
by the University of Michigan Committee for the Care and Use of 
Animals (permit #08194) and were in accordance with the 
National Institutes of Health Guide for the Care and Use of 
Laboratory Animals. All surgery was performed asepticaUy under 
sodium pentobarbital anesthesia, and every effort was made to 
minimize suffering. 

Recordings of rat sleep under baseline conditions and after 24 h 
of REM sleep deprivation wc-re conducted as previously described 
[34]. Briefly, six male Fischer 344 rats (Simonsen Laboratories, 
Gflroy, CA) were used in the study. Under sodium pentobarbital 
anesthesia, the following electrodes were implanted by aseptic 
surgery for the purpose of electroencephalography (EEC) record- 
ing and analysis of behavioral state: 2 superficial cortical electrodes 
(relative to bregma AP: +0.3 and ML: +1.0 for left frontal; AP: -3.0 
and ML: —2.0 for right parietal), 1 deep electrode targeted to the 
left dorsal hippocampus (AP: -3.0, ML: +2.0 and DV: -2.9), 1 
sinus ground and 1 recording wire into each dorsal neck muscle to 
record nuchal electromyography (EMG). After recovery from 
surgery and habituation to the recording apparatus, each rat was 
recorded under the following conditions: (1) 8 h of natural sleep- 
wake behavior, and (2) 8 h of natural sleep-wake behavior, 
preceded by 24 h of REM sleep deprivation. Animals were 
deprived of REM sleep using the multiple platforms-over-water 
method [50]. Recordings commenced 2 hours into the light phase 
of the 12-hour light: 1 2-hour dark cycle (lights on at 6:00am) to 
which the animals were habituated. 

EEG and EMG recordings were analyzed to determine 
behavioral state. Each 8 h recording session was analyzed in 
10 s epochs. Three trained experimenters were blinded to the 
condition of each rat and scored each record for REM sleep, 
NREM sleep and waking [51]. The average pair-wise agreement 
for the three scorers was >0.8 for total sleep, NREM sleep and 
REM sleep. REM sleep was identified by low amplitude and 
desynchronized cortical EEG, synchronized hippocampal EEG in 
the d band (4-9 Hz) and quiescent EMG. Sleep-scored data from 
5 animals were used in this study as the sleep patterning statistics 
(such as mean bout duration, number of bouts and percent time 
spent in each state) of the 6th animal differed by more than 2 
standard deviations from the average sleep patterning statistics of 
the other 5 animals. 

To avoid effects due to recording initiation and circadian 
modulation, sleep-scored data for hours 2-5 of the 8 h recording 
for each rat (n — 5) were used for model development. 

Model formalism 

We constructed the single flip-flop and coupled flip-flop models 
using our previously developed neuronal population firing rate and 
neurotransmitter formalism [24] . In this formalism, the firing rate 
of a pre-synaptic population, fX(t) (in Hz), induces expression of 
neurotransmitter concentration, cX(t), which, in turn, acts as 
input to post-synaptic populations. For a single flip-flop pair of 
populations, the firing rate and neurotransmitter equations are 

fx' ={X^(gYXcY)-fX)lxx, cX' ={cX^{fX)-cX)l%,x (1) 
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fY' = {Y^{gx,ycX)-fY)hY, cY' ={cY^(fY)-cY)/x,Y (2) 

where X,Y= W,S for the wake/sleep flip-flop or X,Y = R?",R'ff 
for the REM-on/REM-o£F flip-flop. Since we consider only 
inhibitory influences of all neurotransmitters, the weighting 
parameters for the post-synaptic influence of all neurotransmitter 
gij are negative. The time constants T;i',Ty,Tc;i',and x^y reflect the 
time dynamics associated with firing rate response and neuro- 
transmitter expression, respectively, at the population level. 
Because different experimental techniques lead to differences in 
absolute reported neurotransmitter concentrations, we normalize 
each neurotransmitter concentration between 0 and 1. The 
functions Xrj-A') represent the steady state acti\ation functions 
for population firing rates and are given by the standard sigmoidal 
form of neuronal firing rate models [52,53] 

Xa, (c) = ^ X„a.{^ + tanh ((c - Px)l^x)), (3) 

where ax and govern the slope and half-activation level of the 
activation function, respectively. The functions cXaa(S) represent 
the steady state expression functions for neurotransmitter concen- 
trations and follow a saturating profile 

cXr^(f)=iiinh(f/yx), (4) 

with slope dictated by yx- In simulations of the model networks, 
states of wake, NREM sleep, and REM sleep were interpreted 
based on firing rates of neuronal populations and concentration 
levels of their associated neurotransmitters. 

In each flip-flop, transitions between states are governed by a 
homeostatic drive that promotes sleep in the wake/sleep flip-flop 
and promotes REM sleep in the REM-on/REM-off flip-flop. The 
sleep-promoting homeostatic drive, h{t), describes the universally 
recognized propensity for sleep that increases during time awake 
and decreases during sleep, and is thought to involve the 
neuromodulator adenosine (reviewed in [29]). The REM sleep- 
promoting homeostatic drive, stp(t), represents the proposed 
short-term process underlying REM sleep homeostasis [15], and 
while its underlying physiological mechanisms have not been 
identified, we model it similarly as the sleep-promoting homeostat 
such that stp increases during NREM sleep and decreases during 
REM sleep. The equations governing the homeostatic drives are: 

li = HifW - 9w}(hmax - h)/ Zh,up 

,down 

Stp = Hidspn -fR°")(stp„,ax - Stp)/ Xslp^up 

(6) 

-I- HifR"" - Ogon ){stpmm -Stp)/ Xstp,down 

where H(z) is the Heaviside function defined as H(z) = 0 if z<0 
and H(z)=l ifz>0, and Ox {X= W,R'"') arefX threshold values 
indicating the occurrence of wake or REM sleep, respectively. The 
parameters x^m^Xmax {x = h,stp) give the minimum and maximum 
values, respectively, that the drive variables can attain and 
Xx,uf,'^x,do\i:n dictate the time scales of their increase and decrease, 
respectively. To incorporate the sleep-promoting homeostatic 
drive into the wake/sleep flip-flops, we model the effects of 
adenosine on the VLPO [54—56] by setting the half-activation 



level of the sleep population's steady-state activation function, iSoo, 
to be dependent on h. Since there is less physiological evidence to 
support how the REM-sleep homeostatic drive affects either the 
REM-on or REM-off populations, we turn to a previous analysis 
of flip-flop models for REM sleep control that indicated that 
transition dynamics ar<; more robust when the REM sleep- 
promoting homeostatic drive acts on the REM-off population 
[27]. Thus, we set the half-activation level of the REM-off steady 
state activation function, R^ff , to be a function of stp. The 
equations for the sleep and REM-off steady state activation 
functions are as follows (compare to Equation (3)): 

So, (c,h) =\s„^(l+ tanh ((c - Ps{h))/«-s)) (7) 

R'i{c,stp)= \Rf,,{\ + i2.nH{c-li^off{stp))/c^^ff)) (8) 

where fi^^iy) = k^iy — k\) for x = S,R°ff and y = h,stp, respectively. 
For the wake/ sleep flip-flop, fc| < 0 such that as h increases during 
wake, decreases to promote activation of the sleep 

population and terminate the wake bout. For the REM-on/ 
REM-off flip-flop, ^^o^ > 0 such that as stp increases during 
NREM sleep, fi^off{stp) also increases to promote inactivation of 
the REM-off population and allow the REM-on population to 
activate. 

A full listing of model equations and parameter values is given in 
Appendix SI and Table SI. 

Sources of variability in tine model 

To incorporate experimentally-documented variability of neu- 
rotransmitter release into the model [57], we multiplicatively 
scaled the steady state release functions for each neurotransmitter, 
cX^^, by a noise factor, ^x(0 whose amplitude randomly varied 
(with uniform distribution and unit mean) at discrete times 
dictated by a Poisson process. Another source of variability 
included in the model represents external excitatory stimuli to the 
wake population resulting in brief wake bouts that fragment sleep 
states as is characteristic of rat sleep patterning. These inputs are 
modeled by the addition of random amplitude, brief excitatory 
inputs occurring at discrete times dictated by a Poisson process, 
5(t), in the argument of the steady-state activation function of the 
wake population, fFoo. 

Interactions between the wake/sleep and REM-on/REM- 
off flip-flops 

As described above, analysis of the dynamics of state transitions 
in the experimental sleep recordings under baseline conditions and 
during recovery sleep following REM sleep deprivation informed 
the interactions between the wake /sleep and REM-on/REM-off 
flip-flops included in the model (Fig. 4A). Activity of the wake 
population influences the REM sleep homeostatic drive stp such 
that when the model is in the wake state, stp decays to a low level, 
stpr, forcing the half-activation of the REM-off population, 
flj^ofj (stp), to a low level and promoting its activation. When the 
wake population is inactive, stp dynamics are governed by Eq. (6) 
(see Eq. S15 in Appendix SI). We model an indirect effect of 
activity in the REM-on population on the wake population such 
that when the model is in REM sleep, the frequency of the 
external excitatory stimuh targeting the wake population, S{t), is 
increased. 
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Numerics 

Statistical analysis was performed in the MATLAB software 
package (The MathWorks, Natick, MA). The administrative code 
was also \\Tittc-n in MATLAB. Differential equation integration 
was performed using the XPP software package [58] (www.math. 
pitt.edu/bard/xpp/xpp.html), interfaced with MATLAB via 
system calls. Two dimensional bifurcation diagrams were created 
using XPP, and the higher dimensional diagram in Fig. SIB was 
created using MATCONT (www.matcont.ugent.be). 

Supporting Information 

Figure SI Modulation of the single wake/sleep flip-flop 
model by neurotransmitter expression variability. Mod- 
ulation of the single wake/ sleep flip-flop model by neurotransmit- 
ter expression variability. A: h coordinates of the right (RKfy) and 
left (LKw) knees of the hysteresis loop as is varied (black dotted 
curve, t^p</=l) and as £,\y is varied (red solid curve, ifs = l). B: 
Distance between RKw and LKw as both ^5 and £,\y are varied. 
Solid colored bands are contours indicating equal distance 
between the knees. White x's indicate values where the knees do 
not exist. 
(EPS) 

Figure S2 Effects of model variability on distributions 
of wake and sleep bout durations. Eflects of model 
variability on distributions of wake (A,C) and sleep (B,D) bout 
durations. A,B: Variable neurotransmitter expression levels result 
in positively skewed, Gaussian-like distributions with a positive tail 
of long wake bouts. C,D: Random excitatory inputs to the wake- 
promoting population result in a bi-modal distribution of wake 
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Figure S3 Modulation of the REM-NREM hysteresis 
loop in the post-REM sleep deprivation condition. Values 
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(lower curve, LKg) of the S-shaped curve of fixed point solutions 
(5) of the fast subsystem of the REM/NREM flip-flop, as a 
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(\crti{ al distance between curves) grows as fc^„^ decreases from the 
baseline value of 7. 
(EPS) 

Appendix SI Full model equations, analysis of the 
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conditions. 
(PDF) 

Table SI Model parameter values for baseline sleep 
case. 

(PDF) 

Autiior Contributions 

Conceived and designed the experiments: JRX) GAM VB. Performed the 
experiments: JRD GAM VB. Analyzed the data; JRD VB. Wrote the 
paper: JRD VB. 



16. Amici R, Cerri M, Ocampo-Garces A, Baracchi F, Dentico D, et al. (2008) Cold 
exposure and sleep in the rat: Rem sleep homeostasis and body size. SLEEP 3 1 : 
708-715. 

17. Saper CB, Chou TC, Scammell TE (2001) The sleep switch: hypothalamic 
control of sleep and wakefulness. Trends Neurosci 24: 726-31. 

18. Saper CB, Scammell TE, Lu J (2005) Hypothalamic regulation of sleep and 
circadian rhythms. Nature 437: 1257-63. 

19. Fuller PM, Saper CB, Lu J (2007) The pontine REM switch: past and present. 
J Physiol 584; 735-741. 

20. Tamakawa Y, Karashima A, Koyama Y, Katayama N, Nakao M (2006) A 
quartet neural system model orchestrating sleep and wakefulness mechanisms. 
J Neurophysiol 95; 2055-69. 

21. Phillips AJ, Robinson PA (2007) A quantitative model of sleep-wake dynamics 
based on the physiology of the brainstem ascending arousal system. J Biol 
Rhythms 22: 167-79. 

22. Diniz Behn CG, Brown EN, Scammell TE, Kopell K] (2007) .\ mathematical 
model of network dynamics governing mouse sleep-wake behavior. 
J Neurophysiol 97: 3828-40. 

23. Rempe MJ, BestJ, Terman D (2010) A mathematical model of the sleep/wake 
cycle. J Matii Biol 60: 615-644. 

24. Diniz Behn C, Booth V (2010) Simulating microinjection experiments in a novel 
model of the rat sleep-wake regulatory network. J Neurophysiol 103; 1937—1953. 

25. Kumar R, Bose A, Mallick BN (2012) A mathematical model towards 
understanding the mechanism of neuronal regulation of wake-nrems-rems 
states. PLoS One 7: e42059. 

26. Robinson PA, Phillips AJK, Fulcher BD, Puckeridge M, Roberts JA (2011) 
Quantitative modelling of sleep dynamics. Philos Trans A Math Phys Eng Sci 
369: 3840-3854. 

27. Diniz Behn C(j, Ananthasubramaniam \, Booth V (2013) Contrasting existence 
and robustness of rem/ non-rem cycling in physiologically based models of rem 
sleep regulatory networks. SIAM J Appl Dyn Syst 12: 279-314. 

28. Lo GC, Chou T, Penzel T, Scammell IE, Strecker RE, et al. (2004) Common 
scale-invariant patterns of sleep-wake transitions across mammalian species. 
PNAS 101: 17545-17548. 

29. Porkka-Heiskanen T (2013) Sleep homeostasis. Curr Opin Neurobiol 23; 799- 
805. 

30. Rinzel J (1985) Bursting oscillations in an excitable membrane model. In: 
Sleeman B, Jarvis R, editors, Ordinary" and Partial Differential Equations; 
Proceedings of the 8th Dundee Conl'erence, Lecture Notes in Mathematics 
1151, Springer, pp. 304-316. 



PLOS ONE I www.plosone.org 



15 



April 2014 | Volume 9 | Issue 4 | e94481 



Coupled Flip-Flop Model for REM Sleep Regulation 



31. Ermcntrout G, Terman D (2010) Mathematical Foundations of Neuroscience. 
Springer. 

32. Izhikevich E (2007) Dynamical Systems in Neuroscience; The Geometry of 
Excitability and Bursting. The MIT Press. 

33. Diniz Behn C, Booth V (2012) A fast-slow analysis of the dynamics of REM 
sleep. SIAM J AppI Dyn Syst 11: 212-242. 

34. Mashour OA, Lipinski WJ, Matlcn LB, Walker AJ, Turner AM, et al. (2010) 
Isourane anesthesia docs not satisfy the homeostatic need for rapid eye 
movement sleep. Anesth Analg 110: 1283-1289. 

35. Bcnington J, HeUer H (1994) REM-sleep timing is controlled homeostatically by 
accumulation of" REM- sleep propensity in non-REM sleep. AmerJ Physiol 266: 
1992-2000. 

36. Benington J, Heller H (1994) Does the hinction of REM sleep concern non- 
REM sleep or waking? Prog Neurobiol 44: 433^49. 

37. Diniz Behn GG, Kopell N, Brown EN, Mochizuki T, Scammell TE (2008) 
Delayed orexin signaling consolidates wakefijlness and sleep: physiology and 
modeling. J Neurophysiol 99: 3090-3103. 

38. Chervin RD, Fetterolf JL, Ruzicka DL, Thelen BJ, Bums JW (2009) Sleep stage 
dynamics differ between children with and without obstructive sleep apnea. 
SLEEP 32: 1325-1332. 

39. Bianchi MT, Cash SS, Mietus J, Peng CK, Thomas R (2010) Obstructive sleep 
apnea alters sleep stage transition dynamics. PLoS One 5: el 1356. 

40. Klerman EB. Wang W, Dully JE, J DD, Czeisler C, et al. (2013) Survival 
analysis indicates that agc-rclatcd decline in sleep continuity occurs exclusively 
during nrem sleep. Neurobiol Aging 34: 309-318. 

41. Doksum KA, Sievers GL (1976) Plotting with confidence: Graphical compar- 
isons of two populations. Biometrika 63: 421^34. 

42. Chu-Shore J, Westover MB, Bianchi MT (2010) Power law versus exponential 
state transition dynamics: application to sleep-wake architecture. PLoS One 5: 
el 4204. 

43. Morris C, Leear H (1981) Voltage oscillations in the barnacle giant muscle fiber. 
BiophysJ 35: 193-213. 

44. Fleshner M, Booth V, Forger DB, Diniz Behn G (20 1 1) Multiple signals from the 
suprachiasmatic nucleus required for eircadian regulation of sleep-wake 
behavior in the nocturnal rat. Phil Trans Roy Soc A 369: 3855—3883. 



45. MontiJM, Torterolo P, Lagos P (2013) Melanin-concentrating hormone control 
of sleep-wake behavior. Sleep Med Rev 17. 

46. Jego S, Glasgow SD, Gutierrez Herrera D, Ekstrand M, Reed SJ, et al. (2013) 
Optogenetic identification of a rapid eye movement sleep modulatory circuit in 
the hypothalamus. Nat Neurosci 16: 1637—1646. 

47. de Andres I, Garson M, Villablanca J (2003) The disconnected brain stem does 
not support rapid eye movement sleep rebound following selective deprivation. 
SLEEP 26: 419-425. 

48. Sakurai T (2007) The neural circuit of orexin (hypocrelin): maintaining sleep 
and wakefulness. Nat Rev Neurosci 8: 1 71-- 181. 

49. Borbcly A (1982) A two process model of sleep regulation. Hum Neurobiol 1: 
195-204. 

50. Jouvet D, Vimont P, Delorme F, Jouvet M (1964) Study of selective deprivation 
of the paradoxical sleep phase in the cat. C R Seances Soc Biol Fil 158: 756—759. 

51. Gross B, Walsh CM, Turakhia AA, Booth V, Mashour GA, et al. (2009) Open- 
source logic-based automated sleep scoring software using electrophysiological 
recordings in rats. J Neurosci Methods 184: 10-18. 

52. Wilson HR, Cowan JD (1972) Excitatory and inhibitory interactions in localized 
populations of model neurons. BiophysJ 12: 1—24. 

53. Dcco G, Jirsa VK, Robinson PA, Breakspear M, Friston K (2008) The dynamic 
brain: from spiking neurons to neural masses and cortical fields. PLoS Comput 
Biol 4: el000092. 

54. Morairty S, Rainnie D, McCarley R, Greene R (2004) Disinhibition of 
ventrolateral preoptic area sleep- acti\'e neurons by adenosine: a new mechanism 
for sleep promotion. Neuroscience 123: 451-7. 

55. Chamberlin NL, Arrigoni E, Chou I'C, Scammell TE, (irecnc RW, et al. (2003) 
Effects of adenosine on GABAergic synaptic inputs to identified ventrolateral 
preoptic neurons. Neuroscience 119: 913-8. 

56. Gallopin T, Luppi PH, Cauli B, Urade Y, Rossier J, et al. (2005) The 
endogenous somnogen adenosine excites a subset of sleep-promoting neurons via 
A2A receptors in the ventrolateral preoptic nucleus. Neuroscience 134: 1377—90. 

57. Aston-Jones G, Bloom FE (1981) Activity of norepinephrine-containing locus 
coeruleus neurons in behaving rats anticipates uctuations in the sleep-waking 
cycle. J Neurosci 1: 876-886. 

58. Ermcntrout G (2002) Simulating, Analyzing, and Animating Dynamical 
Systems: A Guide to XPPAUT for Researchers and Students. SIAM. 



PLOS ONE I www.plosone.org 



16 



April 2014 | Volume 9 | Issue 4 | e94481 



